function CTPipeline(Img, info)
% CT pipeline of Xglom
%   The pipeline takes a 3D image (assuming the image is in the local cache)
%   under its Daris ID.  CT images come in DICOM or TIFF stacks.
%   Cleaning the CT image involves gathering a histogram of image values. 
%   Extrema in the histogram indicate transitions (minima) or concentrations
%   of materials in the 3D image. The CT TIFF stack has empty space, tube,
%   embedding material, un-infused tissue, mildly infused tissue, 
%   strongly infused tissue.
%   
%
%  For testing purposes use:
%  >>  [img info] = GetCTtiffImage('./Cache/1008.2.18.7.1.1.2')
%  >> CTPipeline(img,info)
%
%   - Michael Eager,   (michael.eager@monash.edu)


%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
% 
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.



addpath(fullfile(pwd,'plugins','progressbar'));
handles.debug = 0;
handles.DatasetID='1008.2.18.7.1.1.2';
ShowStatus(['Running Xglom CT Pipeline version 1.0 ... '],handles);
tstart=tic;

if nargin < 2
ShowStatus('Xglom  CT Pipeline version 1.0: Loading Tiff image',handles);
[handles.Img handles.Info ]= GetCTtiffImage(fullfile(pwd,'Cache',handles.DatasetID));
else
    handles.Img = Img;
    handles.Info = info;
end

data_cell_string = fileread(fullfile(pwd,'Cache',handles.DatasetID,'info.txt'));
data_cell_string = regexp(data_cell_string,'\s+','split');

for i=1:length(data_cell_string)
    if ~isempty(strfind(data_cell_string{i},'NZ='))
        handles.NumSlides =  str2double(data_cell_string{i}(4:end));
    end
    if ~isempty(strfind(data_cell_string{i},'DX='))
        handles.PixelResolution =  str2double(data_cell_string{i}(4:end));
    end
    if ~isempty(strfind(data_cell_string{i},'DZ='))
        handles.SliceSeperation =  str2double(data_cell_string{i}(4:end));
    end
    if ~isempty(strfind(data_cell_string{i},'TX='))
        handles.XTiles =  str2double(data_cell_string{i}(4:end));
    end
    if ~isempty(strfind(data_cell_string{i},'TY='))
        handles.YTiles =  str2double(data_cell_string{i}(4:end));
    end
end



nimg=NormaliseImage(single(handles.Img));

midZaxis = floor(size(nimg,3)/2);



gui_active(1);      % will add an abort button
hDummy=figure;
hProgress=figure;
hProgress_0 = hProgress;
hProgress  = progressbar([],0,'CT Pipeline Phase 1 of 2: Image Cleaning Procedure (8 steps)' );
ShowStatus('Xglom CT Pipeline version 1.0: Phase 1 of 2: Image Cleaning Procedure (8 steps)',handles);
incr_progress = 0.01; %double(1)/double(60); %estimate 240 secs for imfill, 50 reps for bwconncomp
%disable original hProgress figure
close(hProgress_0);
close(hDummy)


%% Phase 1a: Grab grayscale mask
ShowStatus('CT Pipeline Phase 1 Step 1 of 6: Measure image statistics',handles);


% calculate a histogram of the image values
midslice=squeeze(nimg(:,:,midZaxis-50:midZaxis+50));

Vreduced = reducevolume(nimg,[10 10 10]);
tic;[hRed xRed] = hist(Vreduced(:),300);toc
hRedsmooth = smooth(hRed,20);

tic;[hOrig xorig] = hist(midslice(:),300);toc
clear midslice;
% smooth the histogram
hsmooth = smooth(hOrig,20);

% find the maxima and minima
out = extrems(hsmooth);

minextrems = out.minx( find( (xorig(out.minx) > 0.05) &   (xorig(out.minx) < 0.5) ));
maxextrems = out.maxx( find( (xorig(out.maxx) > 0.1) &   (xorig(out.maxx) < 0.5) ));

if numel(minextrems) >= 2 && numel(maxextrems) >= 2
    
    
    % find the two minima indicating the tube boundary
    % tube_boundary = x(out.minx( find( (x(out.minx) > 0.25) &   (x(out.minx) < 0.5) )))
    tube_boundary = xorig(minextrems(1:2));
    tube = xorig(maxextrems(1));
    tissue = xorig(maxextrems(2));
    
else
    error('CT pipeline unable to calculate reasonable image ranges in histogram');   
end


telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])


if (handles.debug == 1 )
  % x, log10(h),'-b', 
  hHist=figure; plot(xorig,log10(hsmooth),'-k');
  hold('on')
  plot(xorig(out.minx),log10(out.miny),'*r','Markersize',10)
  plot(xorig(out.maxx),log10(out.maxy),'*g','Markersize',10)
  hold('off')

    fig=figure;imshow(nimg(:,:,midZaxis)); 
    %pause; 
    uiwait(handles.figure1);
    close(fig);close (hHist)
    %   ShowImage(maske,handles);
    %   Continue();
end
if ~hProgress
    return
end

clear hOrig

%% Phase 1b: 
ShowStatus('CT Pipeline Phase 1 Step 2 of 6: Finding embedding solution region. (Est. time 6 minutes.)',handles);
hProgress = progressbar( hProgress, incr_progress*4, 'CT Pipeline Phase 1 Step 2 of 6: Finding embedding solution region. (Est. time 6 minutes.)');


CC = bwconncomp( (nimg>tube_boundary(2)));
numPixels = cellfun(@numel,CC.PixelIdxList);
[xPix idx]=sort(numPixels,'descend');

%Elapsed time is 298.715674 seconds.


telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

if ~hProgress
    return
end

%% Phase 1c: Fit bounding box around region and create cropped image
ShowStatus('CT Pipeline Phase 1 Step 3 of 6: Bounding box of the biggest region',handles);
hProgress = progressbar( hProgress,incr_progress*40, 'CT Pipeline Phase 1 Step 3 of 6: Get the bounding box of the biggest region (i.e. contents of tube). (Est. time 40 sec) ');


% Get the bounding box of the biggest region -> contents of tube
s = regionprops(CC, 'BoundingBox');
BB = s(idx(1)).BoundingBox;

% Elapsed time is 37.830972 seconds.
clear CC s


% crop the image and normalise values
cropped_img = NormaliseImage(single(nimg(floor(BB(1)):ceil(BB(1)+BB(4)),...
    floor(BB(2)):ceil(BB(2)+BB(5)),floor(BB(3)):ceil(BB(3)+BB(6)))));
clear nimg
midZaxis = floor(size(cropped_img,3)/2);

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

%imshow(squeeze(Ilog(midaxis,:,:)))
if handles.debug
  fig=figure;
  imagesc(squeeze(cropped_img(:,:,midZaxis)));
  
  %pause; 
  uiwait(handles.figure1);
  close(fig);
end


%% Phase 1d:  Smooth out Ilog greater than zero with inhomogenityCorrection
ShowStatus('CT Pipeline Phase 1 Step 4 of 6: Crop the CT image.',handles);

hProgress = progressbar( hProgress,incr_progress*3, 'CT Pipeline Phase 1 Step 4 of 6: Crop the CT image.');


crop_range_img = cropped_img;
crop_range_img(cropped_img < tissue) = tissue;
% crop_range_img(cropped_img>0.6) = 0.6;
crop_range_nimg = NormaliseImage(crop_range_img);

clear crop_range_img 
midZaxis = floor(size(crop_range_nimg,3)/2);


%tic;corrMask=inhomogenityCorrection((Ilog.^2).*(Ilog>0),3./60);toc

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

%imshow(squeeze(corrMask(midaxis,:,:)))
if handles.debug ==1
    fig=figure;imagesc(crop_range_nimg(:,:,midZaxis));
    %pause;
    uiwait(handles.figure1);
    close(fig);
end
if ~hProgress
    return
end



%% Phase 1g: Image filling
ShowStatus('CT Pipeline Phase 1 Step 5 of 6: Laplace-of-Gassian filter.',handles);
hProgress = progressbar( hProgress, incr_progress, 'CT Pipeline Phase 1 Step 5 of 6: Laplace-of-Gassian filter. (Est. time 5 minutes.)');

% create image filter of laplacian-of-gaussian operator
hLoG = fspecial3('log',5);

Ilog_tissue= imfilter((crop_range_nimg).*(cropped_img>(tissue)),hLoG,'replicate');

 telapsed = toc(tstart);
 disp(['Processing time (sec): ' num2str(telapsed)])
% 
% %imshow(squeeze(depth(midaxis,:,:)))
 if handles.debug==1
     fig=figure;
     imagesc(squeeze(Ilog_tissue(:,:,midZaxis))); 
     pause;  
     uiwait(handles.figure1);
     close(fig);
 end
 if ~hProgress
     return
 end
 


%% Phase 1h: Normalise the depth image and save processed data
ShowStatus('CT Pipeline  Phase 1  Step 5 of 6: Normalise LoG output and Save the processed image',handles);
hProgress = progressbar( hProgress, incr_progress*19, 'CT Pipeline  Phase 1  Step 6 of 6: Save the processed image');

Ilog_tmp = (Ilog_tissue.*(Ilog_tissue<0).*(cropped_img>(tissue*1.01)));
Ilog_tmp(Ilog_tissue < -0.05)=-0.05;

depth = 1- NormaliseImage(Ilog_tmp);

%depth = crop_range_nimg;  % NormaliseImage(depth);




telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])
if exist(fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat'),'file')
  system(['chmod ug+w ' fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat')]);
end
save(fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat'),'depth','crop_range_nimg', 'Ilog_tissue','tissue', 'tube_boundary','-v7.3');
if exist(fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat'),'file')
  system(['chmod ug+w ' fullfile(pwd,'Cache', handles.DatasetID ,'depth_mask_01.mat')]);
end
clear crop_range_nimg cropped_img Ilog_tissue

if ~hProgress
    return
end

%% Phase 2: Thresholding
ShowStatus('Phase 2: Performing Threshold Statistics (50 steps)',handles);
hProgress = progressbar( hProgress, incr_progress/10, 'Phase 2: Performing Threshold Statistics (50 steps)');

telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])
ThresholdRange= 0.11:0.01:0.5;
stored_hist=zeros(50,1000);
diams_range=1:0.1:100;
stored_diams = zeros(0,141);
    area_max = 1000; %max(areas);
    area_range=(1:area_max);
    
for counter=1:length(ThresholdRange) % Loop over different Thresholds
    %% Thresholding
    Holes=depth>ThresholdRange(counter);
    
    %% Discard single voxels
    % ConArray=conndef(3,'minimal');
    % ConArray(2,2,2)=0;
    % Holes=imdilate(Holes,ConArray).*Holes;
    % Holes=single(Holes);
    
    %% Statistics, Label Glomerular Regions
    CC = bwconncomp(Holes,6);
%    if counter == 20
%        labeledgloms =  labelmatrix(CC);
%    end
    
    % calculate area histogram
    areas = cellfun(@numel,CC.PixelIdxList);

    [area_hist, area_axis]=hist(areas,area_range);
    stored_hist(counter,:) = area_hist;
    
    
    % Total number of detected objects in area_range
    disp(['sum(area_hist(1:100)):  ', num2str(sum(area_hist(5:200)))]);
    
    % Diameter histogram from area data assuming spherical glomeruli
    [diams x]=hist((6/pi*areas(areas>2)).^(1/3),diams_range); diams(end)=0;
    stored_diams(counter,:) = diams;
    
    %disp(['sum(diams(1:100)):  ', num2str(sum(diams(1:100)))]);
    
    
    if ~hProgress
        return
    end
    if counter == length(ThresholdRange)
        ShowStatus('Phase 2: Completed Xglom pipeline. Saving Data.',handles);
        
        hProgress = progressbar( hProgress, incr_progress, ['Phase 2: Completed Xglom pipeline. Saving Data.']);
    else
        ShowStatus(['Phase 2: Threshold Statistics (step ' num2str(counter) ' of 50)'],handles);
        
        hProgress = progressbar( hProgress,incr_progress , ['Phase 2: Threshold Statistics (step ' num2str(counter) ' of 50)']);
    end
end
telapsed = toc(tstart)


ShowStatus('Calculating optimal threshold for detecting gloms ...',handles);

  % Find the best threshold assuming the gloms are around 100 voxels in size
areas_estimate = sum(stored_hist(:,100:150)');
Areas_minmax= extrems(areas_estimate)

% h100= smooth(areas100,3);
% diams_tmp = stored_diams';
% 
% sumdiams=sum(stored_diams');
% diam_minmax = extrems(sumdiams);
% 
% areas_tmp = stored_hist';
% areas_tmp(end,:) = 0;
% areas_tmp(1:20,:) = 0;
% hSmooth = smooth(areas_tmp(5:300,11),20);
% hSmooth = smooth(hSmooth,20);
% sumhist=sum(areas_tmp);
% Areas_minmax = extrems(hSmooth)
if  ~isempty(Areas_minmax.maxx)
   max_i = Areas_minmax.maxx(1);
   max_d = Areas_minmax.maxy(1);
else
    [max_d max_i] = max(areas_estimate);
end
  optimal_threshold = ThresholdRange(max_i);

  % Find the mean and distribution   
  % this process should be done using a fitting method
%   opt_hist = stored_hist(max_i,:)
%   hSmooth = smooth(opt_hist,10,'lowess',10);hSmooth = smooth(hSmooth,10);
%   GlomAreas_minmax = extrems(hSmooth)
% 
% if  ~isempty(GlomAreas_minmax.maxx) && ~isempty(GlomAreas_minmax.minx)
% 
%  
%    GlomOptimum_mean = GlomAreas_minmax.maxx(1);
%    GlomOptimim_sd = GlomAreas_minmax.maxy(1);
% 
%    GlomOptimum_sum = stored_hist(max_i,GlomAreas_minmax.minx(1): 300)
% 
% else
%     [max_d max_i] = max(areas_estimate);
% end

opt_hist = stored_hist(max_i,:)
hSmooth = smooth(opt_hist,10,'lowess',10);
x=1:numel(hSmooth);
fvoxels = fit(x(:),hSmooth(:),'gauss3')
if fvoxels.b2 < 300 && fvoxels.b2 > 40
    GlomVOptimum.mean = fvoxels.b2;
    GlomVOptimum.sd = sqrt((fvoxels.c2)/2);
elseif  fvoxels.b1 < 300 && fvoxels.b1 > 40
    GlomVOptimum.mean = fvoxels.b1;
    GlomVOptimum.sd = sqrt((fvoxels.c1)/2);
elseif  fvoxels.b3 < 300 && fvoxels.b3 > 40
    GlomVOptimum.mean = fvoxels.b3;
    GlomVOptimum.sd = sqrt((fvoxels.c3)/2);
else
    disp('Fitting procedure failed to find glom size range between 300 and 40 voxels')
    GlomVOptimum.mean = 0;
    GlomVOptimum.sd = 0;
end

x=diams_range(1:300);%
opt_diams = stored_diams(max_i,1:300)
hSmooth = smooth(opt_diams,10,'lowess');
fdiams = fit(x(:),hSmooth(:),'gauss3')
if fdiams.b1 < 15 && fdiams.b1 > 4
    GlomDiamOptimum.mean = fdiams.b1;
    GlomDiamOptimum.sd = sqrt((fdiams.c1)/2);
elseif  fdiams.b2 < 15 && fdiams.b2 > 4
    GlomDiamOptimum.mean = fdiams.b2;
    GlomDiamOptimum.sd = sqrt((fdiams.c2)/2);
elseif  fdiams.b3 < 15 && fdiams.b3 > 4
    GlomDiamOptimum.mean = fdiams.b3;
    GlomDiamOptimum.sd = sqrt((fdiams.c3)/2);
else
    disp('Fitting procedure failed to find glom diam size range between 15 and 4 microns')
    GlomDiamOptimum.mean = 0;
    GlomDiamOptimum.sd = 0;
end


fit_output = evalc('disp(fc)');
% if handles.debug==1
%  plot(stored_hist(max_i,:)); xlim ([1 400])
%   hold('on'); %plot(hSmooth,'-r')
%   plot(GlomAreas_minmax.maxx,GlomAreas_minmax.maxy,'*r','Markersize',10)
%   plot(GlomAreas_minmax.minx,GlomAreas_minmax.miny,'*g','Markersize',10)
%   hold('off')
% end

ShowStatus(['Optimal threshold: ' num2str(optimal_threshold) '  Number of glomeruli: ' num2str(max_d)], handles);
hProgress = progressbar( hProgress, incr_progress, ['Xglom CT pipeline. Getting optimum labeled image. Saving Data.']);

   
    % re-reun the connected components routine to get the labelled  
    CC = bwconncomp((depth> optimal_threshold),6);
    labeledgloms =  labelmatrix(CC);

%% Write output data to file - ensure writability by users and group
system(['[ -f ./Cache/' handles.DatasetID '/output_01.mat ] && chmod ug+w ./Cache/' handles.DatasetID '/output_01.mat']);
save(['./Cache/' handles.DatasetID '/output_01.mat'],'stored_hist', 'stored_diams', 'labeledgloms','ThresholdRange','optimal_threshold', 'GlomDiamOptimum','GlomVOptimum')
system(['chmod ug+w ./Cache/' handles.DatasetID '/output_01.mat']);



ShowStatus(['Xglom Pipeline complete for dataset ' handles.DatasetID],handles);
hProgress = progressbar( hProgress, -1 , 'Image Processing Complete. Processing output plots.');



PlotXGlomOutput(stored_hist,stored_diams,handles.DatasetID,...
    handles.PixelResolution,ThresholdRange);
handles.stored_hist = stored_hist;
handles.stored_diams = stored_diams;

telapsed = toc(tstart);
disp(['Xglom CT Pipeline Processing time (sec): ' num2str(telapsed)])


function ShowStatus(str,~)
disp(str)
